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Combining high-temperature expansions with the recursion method and quantum Monte Carlo 
f ■ simulations with the maximum entropy method, we study the dynamics of the spin-1/2 Heisenberg 

ON ' chain at temperatures above and below the coupling J. By comparing the two sets of calculations, 

their relative strengths are assessed. At high temperatures, we find that there is a low-frequency 
peak in the momentum integrated dynamic structure factor, due to diffusive long-wavelength modes. 
, This peak is rapidly suppressed as the temperature is lowered below J. Calculation of the complete 

dynamic structure factor S(k,u) shows how the spectral features associated with the two-spinon 
continuum develop at low temperatures. We extract the nuclear spin-lattice relaxation rate 1/Ti 
qq , from the u> — » limit, and compare with recent experimental results for Sr2CuC>3 and CuGeOa. 

We also discuss the scaling behavior of the dynamic susceptibility, and of the static structure factor 
S(k) and the static susceptibility x{k)- We confirm the asymptotic low-temperature forms S(n) ~ 
[In (T)] 3/2 and xO) ~ r _1 [ln (T)} 1/2 , expected from previous theoretical studies. 

£j 1 PACS: 75.10.Jm, 75.40.Gb, 75.50.Ee, 76.60.-k 

S ' I- INTRODUCTION 

Quantum antifcrromagnets represent an important class of systems in both theoretical and experimental condensed 
matter physics. In recent years, greatly improved precision of neutron scattering and NMR experiments have made 
possible very detailed comparisons with theoretical predictions. A number of new materials have been synthesized 
which appear to bajiear-perfect realizations of the simple spin-half Heisenberg model in various geometries. For 
fNJ , example, Sr2Cu03,tra SrCu203,u and Si^CuC^C^fH comprise structural copper-oxygen units with magnetic properties 
well described by the Heisenberg model on a single ID chain, two coupled chains, and a 2D square lattice, respectively. 
Sr 2 Cu03 is interesting because it appears to be the most perfect ID spin-half Heisenberg system found so far —with 
an exchange J w 2000K and a 3D ordering temperature Tn ~ f\K. Detailed experimental studies of this systemjlcl as 
well as other quasi-lD materials such as KCUF3E1 and CuGeC>3,El have pointed to the need for more accurate theoretical 
studies of the spin dynamics of the 5 = 1/2 Heisenberg chain. Although this model, defined in standard notation by 
£ — ' the Hamiltonian 

H=jJ2 S i' s i+u (J>0), (!) 
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is perhaps the most studied of the basic interacting quantum many-body models, its finite-temperature dynamic 
properties are not fully understood. The low-temperature (T <C J) behavior is controlled by the T = quantum 
critical point (line of critical points to be exact). The powerful machinery of bosonizaijon and conformal field theory 
enables one to make a number of experimentally verifiable predictions in this regimefllj The high temperature regime 
(T J)-, has been studied numerically by short-time or frequency moment expansions combined with the recursion 
methodE The regime of intermediate temperatures T sa J is the most difficult to study theoretically, but is clearly of 
much experimental and theoretical significance, containing the crossover from the diffusive high-temperature behavior 
to the low-temperature regime dominated by elementary-excitations. Here we study the dynamics at intermediate 
temperatures using the high-temperature expansiaaJHTE) technique and a recently developed "stochastic sexies 
expansion" quantum Monte Carlo (QMC) techniqueE30 (an improved variant of the so called Handscomb's methoo£3). 
We have also numerically diagonalized the Hamiltonian for a chain with 16 spins, which although not large enough 
to give reliable results in general, provides for a good test of the other methods in certain-regimes. 

The HTE method has been extensively used to study static properties of spin-models. t3 Here we combine it with 
the continued fraction (or recursion) methods to calculate dynamic properties at finite temperatures. The QMC 
method used here haSj-also previously been applied to both statics and dynamics of Heisenberg models in several 
different geometries £5E£I Accurate results for imaginary-time dependent correlation functions can be obtained down 
to fairly low temperatures. The maximum-entropy (Max-Ent) methocOO is used for analytic continuation to real 
frequencies. This approach has previously been applied to the spin dynamics of the ID Heisenberg model by Deisz, 
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Jarrell and Coxjlaila who focused mainly on the low-temperature dynamic structure factor and the differences between 
half-integer and integer spin. They also discussed at length the accuracy of the Max-Ent method. Here we find that 
for static properties, results obtained using HTE and QMC agree almost perfectly for T/J > 1/8, below which the 
HTE method becomes unreliable. For dynamic properties, the HTE method performs well for T/J > 0.5, and in this 
regime the results agree well with those of QMC and the Max-Ent method. This gives us confidence that the QMC 
and Max-Ent techniques are reliable at lower temperatures as well. 

Our main results are the following: At high temperatures, we find that the dynamic structure factor S(k,uj) grows 
sharply (perhaps diverges) as fc — > and u> — > 0. This diffusive behavior leads to a low-frequency peak in the 
momentum average 

/dk 
-|^(fc)| 2 5(fc,c), (2) 

if the form factor A(k = 0) is non-zero. In an NMR experiment, depending on the A(k) corresponding to a gipm 
material and nucleus under study, this can have large effects on the spin-lattice relaxation rate, which is given byEJ 

±=2S A (u^0). (3) 

J-i 

QMC+Max-Ent results for 1/Ti were previously reported in Ref. [fa. Here we provide results of higher accuracy, 
obtained by calculating the full momentum dependence of S(k,ui — > 0). 

As the temperature is lowered below T rs J/2, the diffusive peak rapidly diminishes in magnitude, and the low- 
frequency spectral weight shifts to k = ir, as expected. ThCj-QMC results for S(k,u>) clearly show the emergence of 
spectral features that can be associated with the well knownniu'S T = two-spinon continuum. 

Comparisons of the momentum^ and frequency-dependent numerical data with scaling theories at low temperatures 
have been presented elsewhereB 2 ] Here we briefly discuss how the scaling in q/T is violated due to logarithmic 
corrections. We provide numerical results for the temperature dependence of the staggered structure factor S(tt) 
and the static staggered susceptibility x( 7r )- At low temperatures the former behaves as [In (1/T)] 3 / 2 , while the latter 
behaves as T _1 [ln (1/T)] 1 / 2 , both expected from theoretical results. 

In Sec. II we discuss the dynamic structure factor and the computational methods used in this study. The results 
are presented in Sec. Ill, and in Sec. IV we discuss and summarize our main conclusions. 



II. BASIC DEFINITIONS AND NUMERICAL TECHNIQUES 



We begin by reviewing some basic definitions of the static and dynamic correlation functions we wish to calculate. 
Both the neutron scattering intensity and the NMR spin-lattice relaxation rate measure the dynamics of the electronic 
spin system through coupling via the operator St ■ Hence, the relevant dynamic correlation function is (SZk(t)S£ (0)) . 
In a spin-rotationally invariant system, which is considered here, this can be evaluated with respect to any quantization 
axis, and in numerical calculations it is most practical to choose the component diagonal in the representation chosen. 
Hence, we study the time-dependent spin-spin correlation function 

S r (t) = (^(i)S 2 (0)), (4) 

were S*(t) denotes the z-component of a spin-1/2 operator at site r at time t, and brackets denote thermodynamic 
averaging at temperature T/J = /3 _1 . We consider only the case of zero average magnetization; (S*) = 0. The 
dynamic structure factor S(k,u>) is the space-time Fourier transform of Eq. (|I|): 

oo 

S(fc,c) = ]T / dte-^-^iSmsm)- (5) 

Apart from kinematic factors, the neutron scattering cross section is directly proportional to S(k,u>). 

NMR (and related techniques such as NQR) can provide accurate results for the low-frequency dynamics, through 
the spin-lattice relaxation rate, given by Eqs. (B l a nd (||). The hyperfine form-factor A(k) can be obtained from 
the Knight shift, and also from impurity effects. E°rE3 Here we will restrict our attention to the important case where 
the nucleus under consideration resides on the sites of the electronic spins, and assume that the real-space hyperfine 
coupling A(r) has a on-site (direct) term A(0) and a nearest neighbor (transferred) term A(l), giving A(k) = A(Q) + 
2A(1) cos (k). The spin-lattice relaxation rate is then given by 
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±- = 2A 2 (0)S r (lo = lu n ), (6) 

where, R = A(1)/A(0), u>n is the resonance frequency, and we define 

S R (u) = (1 + 2i? 2 )5 (w) + iRSxicj) + 2R 2 S 2 (cj), (7) 

where S r (u>) is the real space dynamic spin correlation function at distance r, i.e., the time Fourier transform of 
Eq -®- 

The static structure factor S(k) is the Fourier transform of the equal-time correlation function, and the static 
susceptibility x(^) is given by the Kubo integral 

P 

xik) = J2e ikr J dr(S z r (T)S*(0)), (8) 
r o 

where S*(t) = e T ^ S^e^ T ^ . S(k) and x(k) can be related to the dynamic structure factor through the sum rulesS 

oo 

S(k) = - I dw(l + e-^)S{k, w), (9a) 
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X(k) = - I duu- 1 [l-e- pu )S{k,u). (9b) 
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Below we briefly describe the numerical techniques we use to calculate the dynamic structure factor of the ID 
Heisenberg model. 

A. High Temperature Expansion and the Recursion Method 

The correlation function (^) has a short-time expansion 

Sr (t) = y2M n ^>-, (io) 

n=0 

where the coefficients M n are defined as frequency moments, 

/OO 1 
^ n S r (Lj), (11) 
-oo <™ 

/oo 
dte~ iult S r {t). (12) 
-oo 

An important related function is the spectral density <J>(w) = (1 + e~P u )S(<jj) 1 which is a real and even function 
of frequency and, consequently, its inverse Fourier transform Co(t) (often called the fluctuation function) has an 
expansion in even powers of time only: 

C (*) = X>l) fe M 2fe (13) 

A short-time expansion is of little help if one is interested in the asymptotic long-time behavior of C'o (t) , unles^ some 
kind of analytic ansatz (most often a gaussian one) is made. To this end, let us define the relaxation functior 



of the time Fourier transform 



f 

co(z)= dte- zt C (t)=Y, M ^ z ~ {2k+1) - ( 14 ) 
J ° k=o 
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From Eq. (|l3|), one then has 



coW = / 75 r— » (15) 

7-00 27T Z + 

and $(a/) = 2 lim e ^ +0 i?e[co(e — iw')]. Thus, upon analytic continuation the relaxation function gives the spectral 
density. A useful property of the relaxation function is that it has a continued fraction representation:!] 

c o( z ) = A ° Al • ( 16 ) 
z — , a 2 

To simplify notation we shall in the following write the continued fraction as 

co(z) = A /(z + A x /(z + A 2 /(z + ...))) • (17) 

The first K of the continued fraction coefficients are uniquely determined by the corresponding first K even frequency 
moments (0) through an iterative procedure described in Ref. ||. Of course, we just traded the short-time behavior 
of Co(t) for the large-z behavior of Co(z), which does not bring us any closer to the desired z ~ region. But, as 
described in detail in the book by Viswanath and Mullerjj the relaxation function is uniquely determined by the 
sequence A^, which contains important information about the asymptotic behavior of cq(z). Namely, for the isotropic 
Heisenberg model (|l|), the A-sequence grows with the index k according to a power law, A^ ~ fc A , with 1 < A < 2, 
which fixes the high-frequency behavior of the spectral density to $(w) ~ exp (— |w| 2 / A ). Moreover, oscillations of 
the odd (or even) continued fraction coefficients around the k x curve contain information on the infrared behavior 
of $(w). The simplest function that incorporates both the high- and low-frequency behavior typical of critical spin 
systems has the form 

= I m?r+ ^^ 1 - r exp f- 1 - |2A ) ■ (18) 

The frequency moments of this function are known to be 



: r(|(l + q + 2fe)) 

r(|(i + «)) 



M 2k = uf - (19) 



and the corresponding continued fraction coefficients A k can be calculated from them numerically. Of course, an 
approximation of the spectral density <&(u>) of the system under study by the model density $(w), with parameters 
ujq, a, A determined from a "given" sequence Afc, would be just marginally better than the often used gaussian ansatz. 
Instead, Muller and collaborators (see Ref. [)] and references therein) devised a more accurate procedure, which we 
describe here for completeness. 

Suppose that we have calculated the first K even moments of the true spectral density $>(uj). Then we calculate the 
corresponding sequence A&, and try to approximate it by the model sequence A& by minimizing the sum ^2 k=k (A& — 
A*,) 2 with respect to the parameters woi o:, A. The lower cut-off k m i n (= 3 in our study) is necessary because the 
first few coefficients Aq, Ai, Afc min tend to deviate significantly from the asymptotic behavior represented by 
Having determined the parameters of the fit we may find exactly the A"-th level terminator Tk(z) of the model 
relaxation function corresponding to $(w) by (numerically) inverting the equation 

c Q (z) = A /(z + Ai/(z + ... + A K /(z + T K (z)))). (20) 

The terminator thus incorporates information on the asymptotic behavior of the A^-sequence. 
The relaxation function cq(z) in Eq. ( |l7| ) is then approximated as 

co(z) = A /(z + Ax/(z + ...A K /(z + f *(*)))), (21) 

and thus, in addition to the correct large- z behavior contained in the first few exactly known continued fraction 
coefficients, through the terminator Tk{z), cq(z) also incorporates the correct small- z behavior extracted from the Afe 
sequence. Analytic continuation z — > —ioj 1 then gives us the spectral function $(cj') in the whole range of frequencies 
lu' . For the model spectral function $(w) of type ([i"8|), such analytic continuation is performed by hand and only 
requires numerical integration of well-behaved functions. 
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To study the spin dynamics at finite temperature we have calculated moments of S r (u>) by the HTE technique. 
It is well knew™ that these moments, Mk, can be expressed in terms of a thermal expectation value Qt,a fc-fold 
commutator.BEl High temperature expansion can be developed for these quantities by the cluster method.Eil In fact, 
using the same set of clusters, the expansions for the fc-th moment will be complete to order f3 k , where N depends 
on the size of the largest cluster considered. We have calculated up to N — 22 for all non-zero moments. The 
equal-time correlation function is calculated to order /J 20 . It is more convenient to do the calculations for the scaled 
function co(z)/Mq, i.e. the one defined by the normalized set of moments {1, M2/M0, M2K /Mo}, and obtain the 
needed function cq(z) by simple multiplication at the very end of calculations. 

The behavior of the first 7 continued fraction coefficients of the spin autocorrelation function (Sq (t) z SL(0)) as a 
function of inverse temperature is shown in Fig. 1 . From the relation between M^k and the A& sequencesa we have 
A = Mq, and hence A = 1/4 irrespective of temperature in that case. Another property of the A^'s that make them 
convenient for numerical computations is that they are all numbers of order 1, whereas the corresponding frequency 
moments grow very fast with the index, for example M\4,(T = 00) = 166988876. Continued fraction coefficients 
with higher index k are related to the moments M<ik which are determined by spin correlations at larger distances, 
and hence are more temperature sensitive. The HTE ceases to work for A7 at j3 ~ 0.75 where it starts to change 
rapidly, and for Ag at [3 ~ 1.5. Since we want to reach the lowest possible temperature, we restrict ourselves to first 6 
(k = 0...5) coefficients of the sequence, which permits analysis up to (3 ~ 2. Of course, the shorter the A& sequence the 
more uncertain becomes the determination of parameters luq, ct, A, and in each particular study a try-and-see approach 
should be used to find a compromise between these two conflicting requirements. We found that at j3 < 0.5 results 
obtained with K — 7 and K = 5 do not differ much, and upto (3 = 2 the sequence A0...A5 is stable and reliable. 

In Fig. 1 we also show the fc-dependence of the A& sequence for Sr(uj) (Eq.(Q)) for R = and R = —0.5. The 
latter has a vanishing form factor at k — and thus has no contributions from the diffusive modes. It is evident that 
the former sequence exhibits an odd-even oscillation, suggesting an infrared singularity, but this is absent from the 
latter sequence. This ability to recognize the presence or absence of diffusive modes at such a simple level shows the 
power of the recursion method. 

The temperature dependence of the parameters of the fit, Eq.(|l^), is shown in Fig. ||. Notice the drastic difference 
in a, the power of the infrared singularity for R = and R = — 0.5: in the latter case it is always zero, whereas in 
the former it decreases rapidly with temperature, demonstrating the suppression of the diffusive behavior at low T. 

The overall quality of the described procedure can be estimated by direct term-to-term comparison of the real 
sequence Aj, with the model one A/., given in Fig.|| for k = 3, 4, 5. Deviations of Aj, from A^ are most pronounced 
at high temperatures for the R = case, where the low- frequency infrared divergence is the strongest, and almost 
disappear for T below J. The deviations at high T can be reduced by working with longer A^-sequences, as was 
mentioned above, but this would significantly reduce the range of applicability of our calculations. The absolute value 
of the discrepancy between A's does not exceed 5% in the worst case. 



B. Stochastic Series Expansions and the Maximum Entropy Method 

The stochastic series expansion QMC methotS'0 is an improved variant of the so called Handscomb's technique.!! 
It is based on importance sampling of terms of the power series expansion of exp(— /3H), which for a finite system at 
finite (3 can be carried out to all important orders, without introducing systematic errors. The current formulation of 
the method is described in Ref. O, and has jaijeyiously been applied to both static and dynamic properties of Heisen- 
berg models in several different geometries .EjuJ Here we briefly review how imaginary-time independent correlation 
functions are evaluated using this technique. 

The Hamiltonian for an iV-site periodic chain is first written as 

N 

£ = -!E(^-^)+^P (22) 
6=1 

where the operators and Hij, are defined as 

= 2( 3 - S£Sf +1 ), (23a) 
H 2 ,b = S^S^ +1 + S^S^. (23b) 

The partition function, Z = Tr{e - ^^}, is transformed into a sum suitable for importance sampling techniques by 
expanding e~" H in a power series, and writing the trace as a sum over diagonal matrix elements in the standard basis 
{\a)} = {\Sf,...,S* N )}. This gives 
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z=E££^P(§)>|nM4 < 24 > 



1=1 

where S n denotes a sequence of index pairs defining the operator string Y[?=i HaiMr 

S n = [at, h] [a 2 ,6 2 ] ...[on, b n ], a { € {1, 2}, h € {1, . . . , N b }, (25) 

and ri2 denotes the total number of index pairs (operators) [a,, 6,-] with a$ = 2. For even A" (or any TV for an open 
chain) n% is even for all non-zero terms in Eq. (|24|). All terms are then positive, and can be stochastically sampled 
using standard Monte Carlo techniques in the space of index sequences and basis states. 

The simulation is carried out using two different elementary modifications of the index sequence, taking into 
account the constraints imposed by the fact that the operators defined in Eqs. ( p3[ ) are allowed to operate only on 
antifcrromagnetically aligned spin pairs, and that the operator product corresponding to S n must propagate the state 
| a) onto itself. The power n is changed by inserting or removing single diagonal operators [1,6], and the number of spin 
flipping operators is changed by pair-substitutions [1,6], [1,6] <-> [2,6], [2,6] (the two operators selected for updating 
are typically not adjacent in the sequence). The grand canonical ensemble with fluctuating total magnetization, 
m z = J2i S*r can be studied for T / J > 0.08 by also performing spin flips in the states. At lower temperatures the 
acceptance rate for such flips becomes very low, and it is then more convenient to study the canonical ensemble with 
m z — 0. 

QMC calculations can access the dynamic structure factor only through the corresponding correlation function in 
imaginary time, 

C rur2 (r) = {S* ri (T)S* r2 (0))- (26) 

In the stochastic series expansion method, such a correlation function is estimated by measuring the correlations 
between states |a(p)) = \S* \p], . . ■ , Sfj\p]) obtained by propagating |a) in Eq. (24) with p of the operators in the 
product: 

p 

\a(p))=l[H aiM \a), |a(0)) = |a). (27) 
i=i 

The expression for the imaginary-time dependent spin-spin correlation function i E 

I " T m( a _ \n—m n \ _ \ 

<Wr) = ( £ — ^ Vr^^H . (28) 

\ ' p n (n — m)\m\ / 

x m— v ' ' 

where 

n 

c ru rM) = rrrE^b]^b + 4 (29) 

n p=0 

The periodicity of the propagated states imply that [p + n] = S 2 [p] . Note that any value of r is accessible, in. 
contrast to worldline methods where r must be an integer multiple of the time-slice width used in the simulation. E£l 
The corresponding static susceptibility—^ <T2 (i.e., the real-space version of Eq. (^)), can be directly obtained by 
integrating over r in (128). The result isEJ 



n— 1 n—1 

/ ' s / 

Xri,r 2 



3 rr(E^>])(E^>])+/3%f^)- m 

' p=0 p=0 ' 1 



vn(n + l)V^ nl "A^ r2L ^V ' n 

v ' p=0 p=0 

The relation between C ri ,r 2 ( T ) an d the dynamic structure factor defined in Eq. (0) is 



S k (r) = - I dLjS(k,u)K(w,T), (31) 
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where Sk(r) is the Fourier transform 



G 



Sk(T) = ±Y, e ~ lk{r2 ~ ri)c ^r a (r), (32) 

the kernel is 



AT 

ri,r 2 



K(w,t) = e -™ + e- (/3 - r) ' i ', (33) 



and S(k,—Lu) = e ^ u) S{k 1 uu). The analytic continuation of the QMC data for S'fc(r), i.e., inversion of Eq. (31), is 
carried out using the Max-Ent method, E3B which we very briefly review here for completeness. 

For a given momentum transfer k, S(k, lu) is paremetrized as a sum of (5-functions at frequencies u> n , n = 1, . . . N^, 

N u 

S{k,Lj) = Y^S n 6(LJ-<j n ). (34) 

n=l 

The "classic" variant of the Max-Ent method used here amounts to determining the coefficients S n that minimize the 
quantity 

Q = x 2 /2 - otE. (35) 

Here x 2 is t ne deviation of the imaginary-time function Sk (t) corresponding to a particular set of spectral weights 
{Si, . . . , S n } in Eq. (QLirom the QMC estimate Sk{r) and its statistical error ct(t), which is available for a discrete 
set of times t\ , . ■ . , tn t 

x 2 = f^lS q (n)-S q (n)f/a(nf. (36) 

8=1 

E is the entropy, defined with respect to a "default model" m(u>). Using a default which is constant for u> > and 
satisfies m(— u>) = e~P u m(u>), the entropy is (with S(k,uj) assumed normalized to unity), 

E = -J2S n ln(S n )K{u> n ,0). (37) 

i=l 

The parameter a in Eq. ( |35| ) is determined iteratively using a criterion derived using Bayesian logic, leading to the 
most probable S(k,u>) compatible with the QMC data and its errors, as discussed in detail in Ref. [l?]. Typically, 
on the order of N u = 100 — 200 frequencies are used in Eq. (|34|). The amplitudes S n then form a smooth curve 
representing the frequency dependence of S(k,u>). 



C. Comparisons with Exact Diagonalization 



Accuracies of calculations of dynamic quantities using the HTE and QMC methods are difficult to assess rigorously. 
Comparing results obtained in the two different ways provides for a good test. However, the results will never agree 
completely (In contrast, for static quantities the results agree perfectly in the regime where HTE performs well, as 
discussed in Sec. IITC), and it is therefore useful to check the results against other calculations as well. 

For a small system, all the eigenstates can be obtained exactly by numerically diagonalizing the Hamiltonian. Using 
the translational invariance and the conservation of the z-component of the spin, the Hamiltonian consists of blocks 
corresponding to all the combinations of the magnetization m z and the momentum k. For a 16-site system the largest 
blocks have 810 states, and can easily be diagonalized on a workstation. The next two appropriate sizes, N = 18 and 
N = 20, have largest-block sizes of 2704 and 9252 states, respectively, and could be studied with some more effort. 
Here we consider only N = 16. 

The exact dynamic structure factor of a finite system is a set of (5-functions with positions given by the energies 
of the excited states, and amplitudes given by the squares of the corresponding matrix elements of the operator 
For a 16-site system, the number of contributing (5-functions is very small at T = Q, and it is not possible 
to carry out a meaningful comparison with the other methods. As the temperature is raised, the number of 8- 
functions with significant weights increases rapidly, and a relatively smooth spectrum can be obtained by using 
some reasonable broadening of the individual (5-functions. The results can of course not be expected to completely 
represent the thermodynamic limit, but at temperatures where the correlation length is much smaller than the system 



7 



size meaningful comparisons should be possible. We here compare HTE, QMC, and exact diagonalization results for 
S(k,w) at k = tt/2. 

Fig. |J shows results of all the methods at temperatures T/J — 1.0 and 0.5. We have chosen to represent the 
exact results as histograms with a bin width A w = 0.1J. The nth bin contains the integrated spectral weight in 
the frequency range [(n — l)A u ,nA u ]. On this frequency scale, the results still have structure due to the finite size, 
but nevertheless exhibit over-all shapes that the other results can be compared with. Indeed, the HTE spectra have 
shapes that very well match the histograms. The QMC + Max-Ent results are somewhat broader and have more 
rounded shapes (less asymmetric), as would also be expected, but still represent quite reasonable distributions of 
spectral weight. The QMC results were calculated for a system of TV = 128 sites. Results obtained for the same size, 
N = 16, as the exact diagonalization, look very similar, which also indicates that finite size effects are small at these 
temperatures (for the momentum considered here). The relative statistical errors of the imaginary-time data used 
were in the range 10~ 4 — 10~ 3 , which is typical for all the QMC results discussed in this paper (the absolute error is 
typically in the fifth decimal place of the result). 

The limit oj — ► is of special interest, as it determines the spin-lattice relaxation rate. The results shown in 
Fig. ^ indicate that meaningful results can be obtained for this quantity. The differences between the HTE and QMC 
+ Max-Ent results are typically 10-20%. We cannot rigorously establish which calculation is more accurate in the 
lij — s- limit, but based on the better agreement with the over-all shape between HTE and exact diagonalization in 
the whole frequency range, we expect HTE to be more accurate in the temperature regime where it performs well 
(T/J > 0.5). At lower temperatures only QMC + Max-Ent results are available, since exact diagonalization also does 
not provide accurate results below T/J < 0.5, especially for u> — ► 0. Based on high-temperature comparisons such 
as those presented here, we expect the error of the QMC + Max-Ent calculations at lower temperatures to be of the 
order 10-20%. 



III. RESULTS 

The HTE method is best suited for studying the dynamics of ^-integrated quantities. Apart from the results 
for k = tt/2 presented above, we therefore discuss HTE results mainly for momentum averages and the spin-lattice 
relaxation rate 1 /T\ . With the QMC method we have calculated the imaginary-time correlation functions needed to 
obtain the full S(k,ui) for systems with up to N = 128, down to temperatures T/J — 1/8. Results for 1/Ti, as well 
as the transverse rate I/T^g, were already presented in Ref. |l5|, where 1/Ti was calculated for slightly larger systems 
by analytically continuing weighted imaginary-time correlation averages for space separations r < 2, corresponding 
to the model form factor discussed in Sec. II. Here S(k,uj) is first calculated for all momenta, and the momentum 
averaging is carried out after the analytic continuations. This method, though much more cumbersome and time 
consuming, should be more reliable for studying long-time tails such as those resulting from spin diffusion at high 
temperatures. 

Below, in Sec. III-A, we first consider various aspects of the frequency and momentum dependence of the dynamic 
structure factor. In III-B we discuss the spin-lattice relaxation rate, and recent experimental results for Sr2Cu03 and 
CuGeOa. In III-C we discuss the scaling behavior of the dynamic susceptibility, and how it is affected by logarithmic 
corrections. We present explicit calculations of the logarithmic corrections to the temperature dependence of the 
staggered structure factor and the staggered susceptibility. 

A. The Dynamic Structure Factor 

Results for Sq(u>) (corresponding to a form factor A(k) = 1), obtained using the HTE technique at different 
temperatures are shown in Fig. [|. As discussed earlier, we expect these results to be reliable down to T/J = 0.5. We 
find that at very high temperatures there is strong ui~ a divergence which diminishes as T is decreased, and spectral 
weight is transferred to a broad peak at a; ~ 1.5. At T = 1, a peak at w ~ 1.5 is evident together with an infrared 
peak, which is still strong at this temperature. However, at T = 0.5 the diffusion peak is almost invisible and most 
of the spectral weight is concentrated at higher frequencies. 

Sq(u>) calculated using QMC and Max-Ent analytic continuation is shown in Fig. ^, for temperatures down to 
T/J =1/8. Results for systems with N = 64 and N — 128 are compared, in order to assess effects of the system 
size. The agreement with the HTE results is quite good for T/J = 1.0 and 0.5. The diffusive to — > peak at T = 1.0 
is somewhat higher and narrower in the HTE result, whereas it is somewhat more pronounced in the QMC results 
at T/J = 0.5. At high temperatures, the oj — > peak height grows with the system size in a QMC calculation, 
since the diffusive contributions are cut-off at the momentum k\ — 2ir/N in a finite system. Considering the intrinsic 
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difficulties of numerical analytic continuation, in particular of QMC data, the agreement between the HTE and QMC 
+ Max-Ent results has to be considered quite satisfactory. 

Apart from the u ~ peak at high temperatures, the differences between the N = 64 and N = 128 results are 
likely mainly due to statistical fluctuations in the imaginary-time data, which are amplified in the real frequency 
spectra. The dominant peak at lu/J rs 1.5 is very pronounced at f3 = 8, and does not exhibit much size dependence. 
The position of the peak can be understood on thc.-b.asis of the T = two spinon continuum. The Bcthe ansata 
solution gives an exact expression for the lower edgeJ^J Lu m i n (k) — (tt/2) J sin (fc). The upper bound is known to betH 
w m ax(^) = 7rJsin (fc/2). Since the dominant spectral weight is concentrated at the lower edge, one can expect a 
maximum in the momentum average Sq(lu) at lo — tt/2, arising from momenta q ss tt/2, where the lower bound has 
the smallest g-dependence. This is indeed the position of the maximum seen in Fig. g at = 8. 

A maximum at lower frequency also develops in Sq(u)) for T < 0.25 J. It is due the gradual change from relaxational 
to propagating behavior for modes with momenta 37r/4 < k < tt. The peak sharpens and moves towards uj = as 
the temperature is lowered (this trend is seen also in results for lower temperatures, not gbpwn here). This maximum 
is expected from the scaling form for the dynamic susceptibility first derived by SchulzEHa ~ tanh(w/2T)_£br 

T <C J), and has also been discussed in the context of neutron scattering experiments on Cu(C6D5COO)2-3D20.Ed 

Fig. m shows the momentum dependence of S(k,ui — > 0) at several temperatures. It is clear that the low-frequency 
weight is strongly peaked near k = at T/J = 1, but shifts to fc = tt as the temperature is lowered. The strong 
increase as k — > at high temperatures is clearly not captured completely in a small system, due to the discreteness 
of the momentum. Apart from this long-wavelength cut-off, there is no size dependence between N — 64 and 128 
within the fluctuations of the data. We have not explicitly calculated statistical errors of the results, but one can get 
an impression of their order from the (rather low) degree of non-smoothness of the curves. As discussed in Sec. II-C, 
there may be some systematical errors present as well, due to the unavoidable bias of the Max-Ent method used for 
the analytic continuation. In particular, this may be the case close to both k = and k — tt, in the neighborhood 
of the points where low frequency weight first starts to appear (i.e. at temperatures where there is an intermediate 
momentum regime with vanishing low-frequency weight). In these regimes S(k,uj — > 0) may be over-estimated due to 
broadening effects, i.e., low frequency weight may be seen in Fig. ^ where in fact the actual modes only begin to have 
significant weight slightly above uj = 0. Away from these regimes, we expect systematic errors of at most 10-20%, as 
discussed in Sec. II-C. 

We now discuss the full dynamic structure factor in the temperature regime where there is the most significant 
shift in spectral weight from k « to k ~ tt. We present N — 128 results for T/J = 1.0, 0.5, and 0.25 graphed in 
two different ways. First, in Fig. ||, we show 3D graphs with curves of S(k,Lu) for each individual k. This clearly 
demonstrates how the narrow k sa peak present at T/J = 1 is significantly reduced at T/J = 0.5, where there is also 
a massive build-up of spectral weight at momenta close to k = tt. The maximum at k = tt is not yet very pronounced 
at this temperature, however. At T/J = 0.25 the k = peak has vanished almost completely, and the distribution 
of spectral weight starts to resemble what would be expected from the T = two-spinon continuum. Again, the 
lack of smoothness along momentum cross-sections gives an impression of the considerable amplification by analytic 
continuation of the very small statistical errors in the imaginary-time data. The concentration of the weight between 
the lower and upper bounds of the two-spinon spectrum is seen more clearly in the plots of Fig. ^. Here the intensity 
is represented by shades of gray in the (k, w)-plane, and the T = bounds are also shown. It is clear that there is 
very little weight above the upper bounds even at high temperatures, whereas there is significant weight below the 
lower bound. 

One may still wonder how well the Max-Ent method-captures the true temperature dependence of S(k,uj). In the 
previous low-temperature calculations by Deisz et al.fB considerable weight was observed below the rigorous lower 
bound even at temperatures as low as T/J = 1/24, and the expected concentration of weight at the lower edge was 
not well reproduced. Above we have shown that our high-temperature results agree well with HTE calculations. In 
order to further investigate the broadening effects due to Max-Ent analytic continuation, we have also carried out 
simulations for the system size and temperature considered by Deisz et al. (N — 64, j3 = 24). Fig. [n] shows our 
result for k = 3tt/4, which can be compared with Fig. 9 of Ref. [n]. Our result indeed peaks at the lower bound, and 
is significantly less broadened towards lower frequencies. This probably reflects a higher accuracy in the underlying 
imaginary-time data. The broadening that is present at j3 = 24 may still be partly due to temperature effects, but 
is likely mainly Max-Ent induced. This kind of broadening should be a problem primarily in cases where the lower 
edge is sharp, i.e., at very low temperatures. It will be present to some extent also for temperatures and momenta 
where there is very little low-frequency weight, and, as already discussed above, may then lead to an over-estimation 
of S(k,u> — > 0). The broadening below the lowest bound seen in the results of Figs. || and [| is considerably larger 
than in Fig. nO, and we therefore believe that it is mainly due to real temperature effects, with only minor distortions 
due to Max-Ent bias. 







B. Spin-Lattice Relaxation Rate 



Next we discuss our calculations of the nuclear spin-lattice relaxation rates 1/Ti, for different hyperfine couplings 
parametrized by the ratio R in Eq. (|J). 

The QMC results are obtained by averaging the low-frequency dynamic structure factor shown in Fig. [7] (for 
N = 128). As already mentioned, this method differs from the approach previously used in Ref. [l5[ where the 
averaging was done for the short-distance imaginary-time correlation functions, and the analytic continuation carried 
out subsequently. With exact imaginary-time data, the two approaches would yield identical results, but in the 
presence of statistical errors the Max-Ent method will bias the outcome differently in the two approaches. As an 
illustration of this, results for the on-site dynamic structure factor Sq(l>) obtained using the two different orders of 
averaging and analytic continuation are shown in Fig. |ll|. It is clear that the analytic continuation of the on-site 
function Sq(t) misses much of the diffusive behavior for w — > at the temperatures where these contributions are 
the most important. Based on the previous comparisons with HTE results at high temperatures, we believe that the 
approach of continuing S(k, to) individually for each k before averaging is more accurate. This can also be understood 
on grounds that the frequency dependence of S(k,u>), for a given k, has less structure than the momentum average, 
and is therefore easier to reproduce with a high-entropy curve (which is favored by the Max-Ent method). 

The HTE + recursion method is applied to the short-distance correlations, according to Eq. (Q). The zero frequency 
limit will produce a divergent 1/Ti when infrared singularities are present. However, the very low- frequency, long- 
wavelength limit of our results may not be accurate and should be viewed with some caution as the method is based 
on a short-time expansion, involving only finite spin clusters. We have therefore chosen not to focus on the exact form 
of the to — ► 0, q — ► behavior. For calculating 1/Ti we set the nuclear resonance frequency to 0.01J. In a calculation 
using QMC + Max-Ent, the divergence is cut off due to the finite size, and may also be rounded due to resolution 
effects. In real materials, several energy scales including anisotropy and coupling between chains serve to determine 
the cutoff frequency. 

Results for 1/Ti from the recursion method are shown for several different values of R in Fig. |l2], along with the 
QMC + Max-Ent results for R — and R — —0.5. The agreement between the HTE and QMC calculations is quite 
good ioTr-T/J > 0.5, and gives us confidence in the QMC data at lower temperatures. The main difference from the 
previously QMC + Max-Ent results is that the diffusive contributions for R = at T > 0.5 are stronger, as already 
discussed. At low temperatures the present results are «20 percent lower than the previous results, both for R = 
and R = —0.5. It is likely that the present results at low temperatures are still somewhat too high (likely 10 — 20 
percent), due to the broadening effect of the max-ent method discussed above (see also Ref. for a comparison of 
results for the ratio T 2G /VTTi). 

Sachdes: -calculated 1/Ti using the scaling form of S(k,u) obtained from bosonization and conformal field 
theory.BBEl This gives a temperature independent rate, which however is expected to be modified by logarithmic 
corrections to yield a logarithmic divergence as T — ► Qa This is clearly in qualitative agreement with the above results 
for R = —0.5 (where the diffusive contributions neglected in the theory are filtered out), but the accuracy is not high 
enough to extract a form for the low-tera»erature behavior. We have previously discussed other effects of logarithmic 
correction on the NMR relaxation rates.cHl In particular, it can be noted that the ratio Tig I VTTi is modified from the 
universal temperature-independent form, approaching the universal value with infinite slope as T — * (in a manner 
similar to the bulk susceptibility, as discussed inEa). _ 

We now briefly discuss recent experimental results. Takigawa et al. measured both 1/Ti and 1/Tmj in Sr2Cu03,EI 
for which an exchange J ~ 2000K had been previously deduced from susceptibility measurements!!! The hyperfine 
form factor jCpuld be accurately extracted^ using an impurity effect on the NMR line shape predicted by Eggert 
and Affleck.EH Hence, there are no free parameters, and non-ambiguous comparisons with model calculations can be 
carried out. Takigawa et al. concluded that the agreement with the I/T2G calculation of Ref. |l5| was good. For 
1/Ti the lowest temperature studied in Ref. [l^ corresponds to 300K, which was the highest temperature studied 
experimentally. At this temperature, the previous QMC + Max-Ent result was ~ 40% higher than the experimental 
value. As noted above, our improved method of calculating 1/Ti gives a value w 20% lower than before, and hence 
this discrepancy is now largely reconciled (the remaining differences can likely be explained by uncertainties in the 
experimental values of J and the hyperfine couplings, and by remaining errors in the numerical result). As discussed 
above and in Ref. ^|], also the slight temperature dependence observed for T 2 g/VTT\ can be explained theoretically, 
and is largely due to logarithmic corrections. Hence, all the analytical and numerical results are now in good agreement 
with the experiments, indicating that Sr2Cu03 indeed is an almost perfect realization of the spin-half Heisenberg chain 
with only nearest-neighbor interactions. 

We also note that some contributions from diffusive processes, signaled by a dependence of 1/Ti on lvn (i-e., the 
external field strength) were explicitly observed in Sr 2 Cu03E It would clearly be interesting to study this material also 
at higher temperatures, where our numerical results indicate that a considerably stronger diffusive behavior should 
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be observed. 

Another quasi-lD compound which has been studied very actively recently is CuGeC>3, which undergoes a spin- 
Peierls (SP) transition at T w 14K.H This material has a stronger coupling between the chains than S^CuOs, and 
also is expected to have a non-negligible next-nearest-neighbor interaction J2H \/T\ exhibits a strong, almost linear, 
decrease with decreasing temperature, with a reduction by almost a factor of 5 between T « J and T ~ J/5 (If J2 
indeed is significant, J should here— be considered an effective coupling constant) £3ti3 Itoh et al. argued that the 
transferred hyperfine term is small,E3 and hence R w in our notation. Our results shown in Fig. [l^ then indicate 
a reduction by a factor of 2 in the above. temperature regime, which is significantly different from the experiment. 
On the other hand, Fagot-Revurat et alSIi argued that there is a significant transferred hyperfine coupling, and they 
were able to obtain a reasonable fit to the previous 1/Ti results for R « 0.25.E£l Notice that this value of R leads to a 
strong enhancement of contributions from diffusive, k « 0, modes, and as a result, a stronger temperature dependence 
of 1/Ti (as compared with, e.g., the R = case). Numerical calculations for the NMR rates including J2, inter-chain 
couplings, and dynamic phonons have yet to be carried out, and would clearly be very useful for determining the 
importance of additional interactions beyond the standard Heisenberg chain in a realistic model of CuGeOs. 
fact that the temperature dependence of 1/Ti above the SP transition extrapolates to a positive value at T = 
indicates that the frustrating J2 coupling by itself is not sufficient for opening a gap, but it may clearly contribute to 
stabilizing the dimerized state. 
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C. Scaling Behavior 

In this section we discuss the low-temperature scaling behavior of the staggered susceptibility. At a T — quantum- 
critical point with- z = 1, conventional quantum critical scaling implies that the dynamic staggered susceptibility has 
the scaling formEj 

. , a , cq ui N 
xfea; ) = — — ), (38) 

where a is a non-universal number and <f>(x, y) is a universal complex function of both arguments and q measures 
deviations from the antiferromagnetic wave- vector; q = ir — h. This in turn also implies that near the antifcrromagnetic 
wave-vector, the equal-time structure factor S(q) and the antifcrromagnetic susceptibility, xil) have scaling forms 

S(q)/S(0)=f 1 (cq/T), (39a) 
X(q)/X(Q) = h{cq/T). (39b) 

However, for the Heisenberg chain these scaling relations are violated by logarithmic factors.-which are produced by 
marginally irrelevant operators describing interaction between left- and right-moving currentsE-l The failure of scali 
for x{l) was discussed earlier, and a way to take the logarithmic corrections into account analytically was proposed 
Here, in Fig. [l3], we show the scaling plot for S(q). Substantial systematic deviations from scaling in q/T are apparent 
even at the lowest temperatures accessible to us. We also show results for the quantity S'(q) 7 obtained by subtracting 
the ferromagnetic (uniform) component of the spin-spin correlation function from the numerical data for S r (0) before 
Fourier transforming to momentum space. The uniform term is given by — [T/(2csinh(7rTr/c)] _2 , and contains no 
adjustable parameters.!] S'(q) then contains fluctuations with momenta around the antiferromagnetic wave-number 
7r only, i.e. with q ~ 0. However, S' (q) / S' (0) also does not exhibit scaling in q/T, demonstrating the importance of 
subleading non-universal corrections. Unlike what is the case for the static susceptibility,E3 which is dominated by 
low-frequency fluctuations, these deviations cannot be explained by taking only long-distance logarithmic corrections 
into account. 

Turning now to scaling in w/T, we discuss the scaling amplitudes at q — 0, i.e. the staggered structure factor and 
the staggered susceptibility. Including the leading logarithmic corrections, the sum rules @ give the low-temperature 
forms 

S(q = 0) = D s [ln(T s /T)} 3 / 2 , (40a) 
X (q = 0) = D x T- 1 [\n(T x /T)} 1 ^. (40b) 

In Fig. |l4| these quantities are graphed in such a way that the behavior vs. In ( J/T) should be linear if the above 
asymptotic forms hold. The HTE results shown are from several different differential approximants. The agreement 
between the HTE and the QMC data is quite good down to (3 ~ 6 — 8, where the HTE approximants start to deviate 
from each other and from the QMC results, indicating that the HTE method becomes unreliable. Linear behavior 
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consistent with Eqs. ( ^0| ) is indeed observed in the QMC data for J/T > 4. Fitting lines to the results in this regime 
gives D s = 0.094 ± 0.001, D x = 0.32 ± 0.01, T s = 18.3 ± 0.5, and T x = 5.9 ± 0.2. 

The difference in T s and T x may be due to divergent contributions to 5(0) from short distances, careful treatment 
of which would require a separate cutoff. Another important reason is the ferromagnetic contribution to the spin-spin 
correlation function, already discussed above. This contributes to both 5(0) and x(0). In Fig. |l4]we also show results 
for the quantities 5'(0) and x'(0)i obtained by subtracting the ferromagnetic contributions from 5(0) and x(0). The 
parameters obtained in this case are D' s = 0.106 ± 0.001, D' x = 0.34 ± 0.01, T' s = 5.1 ± 0.2, and T' x = 3.9 ± 0.3. The 
high-energy cutoffs T' obtained in this way are now much closer to each other, and also close to the value To = 4.5 
obtained in Ref. 22 by fitting a scaling form to the real-space correlation function at (3 = 32. Within the accuracy of 
the numerical procedures (which apart from statistical errors also include effects from non- asymptotic contributions) 
they can be considered equal. We note that these parameters still should be viewed as 'effective' or temperature- 
dependent, and their .asymptotic zero temperature-values may only be reached at T/J < 0.01 as is the case for the 
uniform susceptibility^.] and the correlation lengthLJ. 

The temperature dependence of x{l = 0) is important in the context of systems of weakhz, coupled spin chains, 
where it determines the critical temperature for ordering into a three-dimensional Neel stateo A recent calculation 
for KCUF3 by Sc1iu1z,e3 using our numerical data for x(l = 0), is in good agreement with experimental results. 



IV. CONCLUSIONS 



In this paper, we have used several numerical methods to study the dynamics of the spin- 1/2 Heisenberg chain at 
intermediate temperatures (above and below J). In order to obtain the dynamic structure factor, we have combined 
high temperature expansions for the frequency moments with the recursion method, and quantum Monte Carlo 
simulations with the maximum entropy method. In some cases, the results of these methods have been compared 
with exact results for a system with 16 spins, as an additional check. 

We find that at high temperatures the HTE + recursion method works very well even for dynamic quantities. Using 
the first six continued fraction coefficients for the relaxation function, we have obtained dynamic susceptibilities at 
various temperatures down to T/J = 0.5. The calculated frequency-dependence of the structure factor at k = ir/2 
agrees remarkably well with the exact results on finite systems. This method is also sensitive to infrared singularities. 
Already at the level of the continued fractions themselves, that is before any numerical extrapolation is done, the 
presence or absence of singularities can be detected. We found that the behavior of the continued fractions for the 
local (fc-integrated) susceptibilities changes qualitatively depending on whether the hyperfine couplings are vanishing 
or non-vanishing at k = 0. However, when infrared singularities are present, their exact forms may not be fully 
captured by these methods. Since the diffusion-related singularities are most robust at infinite T, where they have 
been investigated in the past with higher number of momcnts£j we did not focus on this issue here. In our study, 
the HTE with the recursion method became unreliable below T/J = 0.5 because the extrapolations for the higher 
moments became unstable. By extending the series for the frequency moments, it may be possible to reach still lower 
temperatures also for dynamic quantities. 

The QMC results for the static quantities agree perfectly with the high temperature expansions down to T/J w 
0.1 — 0.2, below which the HTE results become unreliable. Results obtained for the dynamic quantities with the QMC 
+ Max-Ent method are in satisfactory (generally within 10-20 percent) agreement with HTE results above T/J = 0.5 
and, we believe, should have similar reliability down to much lower temperatures. For calculating local quantities 
such as 1/Ti, and incorporating as much as possible the effects of infrared singularities from certain k- regions, it 
appears to be better to carry out the analytic continuation for all individual momenta separately before performing 
the momentum sum, rather than carrying out the analytic continuation for momentum-integrated quantities. 

We have also presented several new results for the Heisenberg chain. These include quantitative estimates for the 
logarithmic temperature dependence of the static staggered susceptibility and structure factor, improved estimates 
for the temperature dependence of the spin-lattice relaxation rate with various choices of hyperfine couplings, as 
well as the full dynamic structure factor at intermediate temperatures. Our results clearly show the shift in the low 
frequency spectral weight from the diffusive modes near k = at high temperatures to the antifcrromagnetic modes 
near k = tx at low temperatures. At low temperatures we observed the development of spectral features that can be 
associated with the two-spinon continuum. Overall, our results at the lowest temperatures are in good agreement 
with with general expectation of these quantities from various analytical studies. Our results also allow us to reconcile 
the measurements of the spin-lattice relaxation rates in Sr2Cu03 with theoretical and numerical calculations for the 
Heisenberg model. _ 

After completing this work, we became aware of a recent paper by Fabricius, Low and Stolze,ciJ discussing exact 
diagonalization results for S(k,u>) of chains with up to 16 sites. We have here used such short-chain calculations 
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mainly as a test of the HTE and QMC + Max-Ent methods in some regimes. Fabricius et alc^ found no signs of 
spin diffusion. Our results indicate that 16 sites is not large enough in general, in particular not for addressing the 
intricate problem of spin diffusion. As can be inferred from Fig. ^, the momentum cut-off ki — ir/8 for a 16-site chain 
prohibits access to most of the long-wavelength regime where,we see signs of diffusive behavior. In the limit T — > 
there is most likely no spin diffusion, as discussed by Sachdev ,H and also supported by the numerical results of Ref. [l5| 
and our present calculations. We stress again that our study is also not accurate enough to resolve the exact form 
of S(k,u) in the k — + 0,lu ~ * limit. We have strong evidence of a sharp peak at high temperatures, but cannot 
determine whether or not it is truly divergent (i.e., whether the long-time behavior of the spin-spin correlations is of 
the standard ID diffusive form ~ t~ x / 2 , or, perhaps, is anomalous). In view of the absence of diffusion in the limit 
T — > 0, a sharp u> — > peak in Sq(uj) at high T, diverging as T — > oo, is perhaps the most likely scenario. In any case, 
as we have discussed above, a sharp maximum should have detectable effects on, e.g., the spin lattice relaxation rate 
in real materials. 
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FIG. 1. Continued fraction coefficients, At, for Sr(u>) as a function of /3 for R = (lower panel) and as a function of the 
index k for R = and —0.5 at /3 = (upper panel). 

FIG. 2. Temperature dependence of the parameters of the fit for two values of R. 

FIG. 3. Comparison of the real and model A-sequences for two values of R at different temperatures. Each panel shows A3 
(bottom), A4 (middle), and A5 (top). 

FIG. 4. Comparison of exact TV = 16 (histograms), HTE (solid curves) and N=128 QMC + Max-Ent (dashed curves) results 
for the dynamic structure factor at k = n/2. 

FIG. 5. The dynamic structure factor Sr(uj) with R = calculated using the HTE method at various temperatures. 

FIG. 6. The dynamic structure factor averaged with a constant form factor calculated using QMC and Max-Ent analytic 
continuation at different temperatures. Dashed and solid lines are results for system sizes TV = 64 and 128, respectively. 

FIG. 7. The low-frequency limit of the dynamic structure factor vs. the momentum. Solid and open circles are for TV = 64 
and 128, respectively. 

FIG. 8. QMC + Max-Ent results for the full dynamic structure factor S(k,u>) at three different temperatures for TV = 128. 
The maxima of the vertical scales are 3.01 (T/J = 1.0), 1.49 (T/J = 0.5), and 1.86 (T/J = 0.25). 

FIG. 9. QMC + Max-Ent results for the dynamic structure factor S(k,ui) at T/J = 1.0, 0.5, and 0.25, with shades of grey 
representing the intensity in the (ui, q) plane. The curves indicate the lower and upper bounds at T = 0. Note that for T = 1.0 
and 0.5, S(k, u>) is sharply peaked at k — > 0, uj — > 0. 

FIG. 10. QMC + Max-Ent results for the dynamic structure factor at k — 37r/4, calculated for a 64-site system at inverse 
temperature f3 = 24. The dashed line is the T = Bethe ansatz lower edge. 

FIG. 11. Comparisons of So(u) calculated by averaging before (dashed curves) and after (solid curves) analytic continuation. 
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FIG. 12. Results for NMR relaxation rates 1/Ti for various values of the hyperfine parameter R, calculated using HTE and 
the QMC and Max-Ent methods. 



FIG. 13. The equal-time structure factor S(q) obtained from HTE and QMC, graphed so that the data would collapse onto 
a single curve if scaling in q/T holds. The QMC results for ji < 8 were calculated in the grand canonical ensemble for N = 256, 
and those for (5 > 16 in the canonical ensemble for N = 1024. 

FIG. 14. The T-dependence of S(k = n) and x(k = n) obtained from HTE and QMC, graphed so that the predicted forms 
( |io| ) give linear behavior. The dashed curves are several different HTE approximant, and the solid circles with error bars are 
the QMC results. The solid lines are fits to the QMC data for T/J < 0.25. The open circles and dotted lines are results and 
fits after subtraction of the ferromagnetic contribution to the spin-spin correlation function. 
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